#include <iostream>
#include "nnec/nnec_utils.h"

using namespace std;
using namespace NNEC;

int main()
{
    cout.precision(10);

    complex<double> a[] =
    {{-5.033858164,486.2389895}, {-9.821994368,-309.0191223}, {-196.405933,-35.09067658},
     {-10.0939426,589.6293652}, {-5.2965558,-5.223079321}, {-5.2965558,-196.405933},
    {-35.09067658,-5.223079321}, {-9.821994368,-309.0191223}, {-5.033858164,486.2389895}};
    complex<double> b[3]={1.0, 1.0, 1.0};
    int p[3]={};
    int d[3]={};

    auto *pa = new complex<double>[9];
    for(int i = 0; i < 9; i++) {
        pa[i] = a[i];
    }

    lu_decomposition<complex<double>>(a, 3, p, d);
    lu_backward_substitution<complex<double>>(a, 3, p, b);

    for(int i = 0; i < 3; i++) {
        for(int j = 0; j < 3; j++) {
            cout << a[i * 3 + j];
        }
        cout << endl;
    }

    cout << endl;

    auto *L = new complex<double>[9];
    auto *U = new complex<double>[9];

    for(int i = 0; i < 3; i++) {
        for(int j = 0; j < 3; j++) {
            if(i == j) {
                L[3 * i + j] = 1.0;
                U[3 * i + j] = a[3 * i + j];
            }
            else if(i > j) {
                U[3 * i + j] = 0.0;
                L[3 * i + j] = a[3 * i + j];
            }
            else {
                L[3 * i + j] = 0.0;
                U[3 * i + j] = a[3 * i + j];
            }
        }
    }

    for(int i = 0; i < 3; i++) {
        for(int j = 0; j < 3; j++) {
            cout << U[3 * i + j];
        }
        cout << endl;
    }

    cout << endl;

    for(int i = 0; i < 3; i++) {
        for(int j = 0; j < 3; j++) {
            cout << L[3 * i + j];
        }
        cout << endl;
    }

    cout << endl;

    for(int i = 0; i < 3; i++) {
        for(int j = 0; j < 3; j++) {
            complex<double> tmp{0.0, 0.0};
            for(int k = 0; k < 3; k++) {
                tmp = tmp + L[3 * i + k] * U[3 * k + j];
            }
            cout << tmp << " ";
        }
        cout << endl;
    }

    cout << endl;

    for(int i : p) {
        cout << i << " ";
    }
    cout << endl;

    cout << endl;

    for(auto & i : b) {
        cout << i << " ";
    }
    cout << endl;

    cout << endl;

    for(int i = 0; i < 3; i++) {
        complex<double> tmp{0.0, 0.0};
        for(int j = 0; j < 3; j++) {
            tmp = tmp + pa[3 * i + j] * b[j];
        }
        cout << tmp << " ";
    }
    cout << endl;

    return 0;
}
